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Many numerical techniques for the description of quantum systems that are coupled to a con¬ 
tinuous bath require the discretization of the latter. To this end, a wealth of methods has been 
developed in the literature, which we classify as (i) direct discretization, (ii) orthogonal polynomial, 
and (iii) numerical optimization strategies. We recapitulate strategies (i) and (ii) to clarify their 
relation. For quadratic Hamiltonians, we show that (ii) is the best strategy in the sense that it 
gives the numerically exact time evolution up to a maximum time fmax, for which we give a simple 
expression. For non-quadratic Hamiltonians, we show that no such best strategy exists. We present 
numerical examples relevant to open quantum systems and strongly correlated systems, as treated 
by dynamical mean-field theory (DMFT). 


I. INTRODUCTION 

Quantum systems coupled to a continuous bath ap¬ 
pear in different fields of physics, such as open quantum 
systems (OQS), strongly correlated many-body physics, 
and spectroscopy and scattering problems. In the con¬ 
text of OQS [1, 2], for instance, a quantum system like 
an atom or a quantum dot is linearly coupled to a contin¬ 
uous bath like a phononic, electronic or photonic reser¬ 
voir, which produces dissipation and decoherence in the 
system. In the context of strongly correlated many-body 
physics, the Anderson impurity model [3] and its general¬ 
izations, which describe clusters of electronic impurities 
coupled to a continuous conduction band of electrons, 
are an important field of study. In addition, they are 
the basis for dynamical mean-field theory (DMFT) [4-6], 
which is the most widely used numerical method to de¬ 
scribe strongly correlated systems in dimensions higher 
than one in physics [7, 8] and is popular also in quantum 
chemistry [9] . A discrete system coupled to a continuum 
appears also in spectroscopy or scattering problems [10], 
leading to a resonance or state with a complex energy 
that due to the imaginary energy component decays in 
time. 

The dynamics of a system that is strongly coupled to 
a continuous environment cannot be described using an¬ 
alytic weak-coupling approaches [1, 2], and requires the 
use of numerical techniques such as exact diagonalization 
(ED), the density matrix renormalization group (DMRG) 
and the numerical renormalization group (NRG). How¬ 
ever, all of these numerical techniques are restricted to 
treating discrete Hamiltonians, and cannot directly deal 
with a Hamiltonian that involves a continuous bath. 
Therefore, it is necessary to construct a discrete approx¬ 
imation to the continuous Hamiltonian. 

In this paper, we analyze the problem of constructing 
the discrete Hamiltonian that best approximates the time 
evolution produced by the continuous Hamiltonian with 
the smallest possible number, W, of discrete degrees of 
freedom. As the many-body Hilbert space grows expo¬ 
nentially with N}), this question is highly relevant, and its 
solution would allow to tackle systems with a complexity 


that is otherwise out of reach. We will show that this 
problem can only be solved for quadratic Hamiltonians. 
For non-quadratic Hamiltonians, we show that no best 
discrete approximation exists, and instead, heuristic ar¬ 
guments have to be used to construct an approximation, 
as already found frequently in the literature [11-18]. 

Let us consider a general setup consisting of a system 
with Hamiltonian iLgys expressed in terms of system op¬ 
erators d) and d (e.g. in the quadratic case iLsys = cod^d), 
which is linearly coupled to a continuous harmonic oscil¬ 
lator bath characterized by a Hamiltonian iLbath, 

H — Hsys T TLbath T LIcoupl; (1^) 

l^b 

Hhe.th= dxxala,,, (lb) 

J a 

fb 

-^coupi — / dx cix h.c., 

J a 

via a “coupling function” V{x). Here, aj, (ax) create 
(annihilate) an occupation of a bath level with energy x. 
This defines the bath spectral density J{x) as [1, 19] 

fb 

J{x) = / dx'\V{x')\^5{x — x') = \V{x)\^. (2) 

J a 

This spectral density, which depends on the continuous 
bath variable x, fully characterizes the influence of the 
bath on the system. Similarly, a system linearly coupled 
to a discrete harmonic oscillator bath is characterized by 
a Hamiltonian 


^discr 

_ TT _|_ rrdiscr j_ rrdiscr 

-^sys “r -^bath -^coupl 
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udiscr _ 

-^bath 

n—1 

(3b) 
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The bath spectral density is a comb of delta peaks and 
not a continuous function as in equation (2) [19], 

Nb 

J^-^--’^(x)=J2\ynfdix-Xn). (4) 

n—1 
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For Nb —>■ oo one can find an iJ^iscr jg equivalent to 
H [20, 21], For Nb < oo, the discrete Hamiltonian (3) can 
only serve as an approximation of the continuous Hamil¬ 
tonian (la). We classify the strategies for constructing 
such an approximation as follows. 

(i) Direct discretization, in which bath energies a;„ and 
couplings Vn are obtained by a discretization of the 
integration interval [a, 5] in (la). This technique is 
standard in the context of NRG [21] and frequently 
used in the context of DMRG [11-18]. 

(ii) Orthogonal polynomials [22], with which the bath 
energies Xn are obtained as the zeros of a polyno¬ 
mial that is associated with a quadrature rule for 
the integration over the continuous bath energies 
X. This has been used in different contexts from 
DMRG to quantum chemistry [20, 23-30]. 

(hi) Numerical optimization, which consists in choosing 
the parameters Xr, and 14? by minimizing a cost 
function [31-33]. 

As strategy (iii) cannot be used to discretize the spec¬ 
tral representation of a bath (see Appendix A), we re¬ 
strict ourselves to strategies (i) and (ii), which we reca¬ 
pitulate in Sec. (HA) and Sec. (HB), respectively. In 
Section (II), we clarify the relation of strategies (i) and 
(ii), which has hitherto been missing from the literature. 
In Sec. (HI), we show that strategy (ii) best describes the 
time-evolution for quadratic Hamiltonians, and that for 
non-quadratic Hamiltonians, there is no such best strat¬ 
egy. Section (IV) presents numerical examples and in 
Sec. (V) we draw the main conclusions of the paper. 


II. RELATION OF DIFFERENT 
DISCRETIZATION STRATEGIES 

Let us introduce the analytic continuation of the bath 
spectral density (2) to the complex plane, the hybridiza¬ 
tion function [19] (see Appendix B) 

A(z) = [ dx — zGC (5) 

Ja Z-X 

with J{x) = jV(a;)j^. By the Sokhotski-Plemelj theorem 
this implies 

J{x) = ——ImA(a; -|- iO). (6) 

TT 

The hybridization function does not contain more in¬ 
formation than J{x) since its real and imaginary parts 
are related by the Kramers-Kronig relation, Re[A(x)] = 
f (ja;= _i r Using the discrete bath 

spectral density J'^““(a;) of (4) to evaluate (5), one ob¬ 
tains 

Adiscr(^) (7a) 

^Z-Xn 

n—1 


A. Direct discretization strategies 


Let us consider the approach of Ref. 25 and rephrase 
the problem of discretizing the Hamiltonian as that of 
discretizing the integral in (5) . The simplest approxima¬ 
tion for an integral is obtained by using a trapezoidal 
integration rule 


A(z) 



\yixn)\'^^Xn _ ^discrj-^-j 

^ Z-Xn ~ 

( 8 ) 


where Xn are linearly spaced node points with spac¬ 
ing Axn- Using this rule to generate an approximation 
Adi^cr^^)^ i.e. demanding the last equality of the preced¬ 
ing equation to hold, it is possible to identify the cou¬ 
plings as 

1V„]2 = \V{Xn)\^NXn (9) 

and the node points Xn as bath energies of (3). 

The strategy using the trapezoidal rule can be im¬ 
proved as follows. Instead of generating a discrete weight 
simply by multiplying the function jU(a:„)j^ with 
the width of the associated interval Nx as in (9), com¬ 
pute the weight 1I4P as an integral of lU(a;)l^ over an 
interval and the bath energies Xn as weighted aver¬ 
ages over this interval 


\Vr,\^ = dx\V{x)\\ (10a) 

dxx\v{x)\'^. (10b) 

This requires to define intervals J„ C [a.6], n= \,...,Nb, 
with InC\ Im = 0 for n ^ m and [a, 6] C 
a linear discretization this generates intervals of equal 
width as in the trapezoidal rule (8). But in general, 
the intervals can have arbitrary widths, and one can 
e.g. dehne a logarithmic discretization, for which the in¬ 
terval widths decrease exponentially for jccj —> 0. This 
guarantees energy scale separation, which is required for 
NRG [21]. ED and DMRG, by contrast, allow for any 
discretization. Within DMRG, for instance, aside from 
the linear [12, 14, 15] and logarithmic discretizations 
[11, 12], it is possible to consider combinations of both 
discretizations [16], combinations of different logarithmic 
discretizations [13], or a cosine-spaced discretization [17]. 
Also, a parabolic discretization has been proposed [18]. 

Within the direct discretization strategy, the discrete 
bath operators in (3) are interpreted as averages of 
the continuous bath operators aj, in (la) over the energy 
interval /„ 

c]) = dxV{x)al. (11) 

The map of i—>■ retains the (anti-)commutation rela¬ 
tion of the continuous operators [a^, a\,]± = 5{x — x') as 
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discretization intervals do not overlap and are normalized 


Limits of the direct discretization strategy 


[Cn,c|„]± — J 

In the context of direct discretization strategies, we 
point out that the discrete representation (3) is typi¬ 
cally referred to as the star representation of the dis¬ 
crete Hamiltonian. This representation is, via a standard 
mapping [ 21 ], unitarily equivalent to a one dimensional 
tight binding chain, i.e. a chain representation (see Ap¬ 
pendix C). This mapping is valid independently of the 
discretization strategy and can even be formally defined 
to map the continuous star Hamiltonian into a chain with 
infinite length [26]. This issue will be further discussed 
in Sect. (HB3). Finally, we note that in the chain rep¬ 
resentation, the logarithmic discretization leads to next- 
neighbour couplings that decay exponentially with the 
distance to the impurity. 


J dx'V*{x)V{x')[a^,al,]± = S„ 


1. New proposals 

To improve the accuracy of the discretization of previ¬ 
ous strategies [11-18, 21 ], it seems reasonable to consider 
a node distribution that uses more nodes in regions where 
the bath spectral weight is larger. Based on this heuris¬ 
tic argument, we propose two different variants of direct 
discretization strategies. 

In the first one, which we refer here simply as the mean 
method, we compute the first bath energy as an average 
over the full support of J{x) 

1 

= ITT —N / dxxJ{x), 

IHotr Ja 

|^"totP= [ dxJ{x). (12) 

J a 

In the next step, we compute X 2 as an average over the 
interval [a,xi], and X 3 as an average over the interval 
[xi, b]. The following steps are repeated in a similar way 
until obtaining energies. Finally, the weights |H„|^ 
are obtained as integrals 

Ax„+x,,+ i)/2 

= / dxJ{x), (13) 

J (Xn-1+Xn)l2 

where for the first (n = 1 ) and the last (n = Nif) integral, 
we replace the lower limit by a, and the upper limit by 
b, respectively. 

Similarly, we define the equal weight method. Here, 
in the first step we define a weight per bath energy 
J^dxJ{x). Then, we define the first interval Ii = 
[a, ai] via 

1 

dxJ{x) = — / dxJ{x), (14) 

J a 



The direct discretization strategies considered in this 
section are based on producing non-equally spaced dis¬ 
cretization intervals to minimize the error of the approxi¬ 
mation A(z) ~ A'^*®'^''(z) for certain values oi z = x + iQ'^, 
i.e. for certain values of the bath energy x. 

The logarithmic discretization, e.g., minimizes the er¬ 
ror in the low-energy limit |x| —>■ 0. This discretiza¬ 
tion then forms a quasi-continuum in a neighborhood of 
X = 0 , and therefore the discretized version of the hy¬ 
bridization in such region is a numerically exact approx¬ 
imation to the continuous one. However, such a good 
approximation for low energies comes at the price that 
for higher energies the discretization becomes crude, and 
the logarithmic approximation is therefore not appropri¬ 
ate to describe the time evolution of the system at short 
and intermediate time scales. Thus, NRG, which uses 
a logarithmic discretization, allows to describe the low- 
energy physics of a system numerically exactly, but gives 
a very rough approximation of high-energy excitations of 
the bath. The proposals described in Sec. (HA 1 ), on the 
other hand, provide a good approximation in those en¬ 
ergy regions where the spectral density is larger in mag¬ 
nitude, which may not necessarily coincide with low en¬ 
ergies. 

In general, none of the direct discretization strategies 
reliably describes the system at all energy scales. More 
precisely, a safe use of these strategies (i.e. unbiased with 
respect to energy) to describe time evolution at short and 
intermediate times scales, requires to consider a relatively 
high number of bath sites {Nb = 30 up to 200, depending 
on the problem [11, 13-18]). 


B. Orthogonal polynomial strategy 

In order to construct a discrete representation of the in¬ 
tegral (5), which is valid for all bath energies x in [a, b], it 
is necessary to use a discretization method in which each 
discretized energy value x„ is computed with information 
of the integrand (5) over the whole integration support 
\a,h]. As will be described in the following, this can be 
achieved by using Gauss-Christoffel type of quadrature 
rules to represent the integral (5), which to our knowl¬ 
edge has for the first time been proposed in Ref. 22. 


1. Gaussian quadrature 


Let us re-express the z-dependent integral (5) in terms 
of the product of a weight function w{x) {w{x) > 0 ) and 
a function f{x,z) (see Ref. 34 for an excellent review on 
the subject). 


and the corresponding first bath energy and weight is 
computed as in (10). The rest of parameters x„ and 
are obtained analogously. 


A(z) = [ dx = [ dxw{x)f{x,z). 

Ja Z-X Ja 


(15) 


z — X 
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Now consider a polynomial interpolant z) of f{x, z) 

with degree iV — 1 (here and in the following, the degree is 
with respect to the argument x, which is the integration 
variable), which is unique and matches f(x, z) at N node 
points Xn, 

f{x,z) = fN{x,z) -\-rN{x,z), (16) 

N 

f n{x^ Z^ = ^ ^ f {Xfi , ~ ^nmt 

n—1 


where ln{x) can be defined as the {N — l)-th order poly¬ 
nomial ln{x) = “ Xm) and 

r^^XjZ) is a remainder. Clearly, if the degree of f(x,z) 
is — 1, one can achieve r^ix, z) = 0 if choosing the N 
node points Xn intelligently, and 


A(z) 

Wn 


.b N 

/ dx W{x)f{x, z) = WnfjXn, z) + Rn{z), 

n=l 


dx w{x)ln{x), 


(17) 


is an exact representation of the integral, i.e. Rn{z) = 0 . 
We refer to Wn as Christoffel weights. It can be shown 
that Rjv(z) = 0 holds even if f{x, z) has a degree smaller 
or equal than 2N — 1, although then rj^{x, z) ^ 0. The 
integration rule is then of degree of exactness 2N — 1. 
The higher the degree of exactness, the smaller is the 
error term i?Ar(z) for the function f{x,z), even if the 
latter has degree higher than 2N — 1. 

To obtain the highest possible degree of exactness 
2N — 1, Posse and Christoffel showed in 1877 that the 
previously referred intelligent choice of the nodes Xn is 
to consider them as the roots of the monic polynomial 
Pn{x) of degree N that pertains to the family of orthog¬ 
onal polynomials obeying 


b 

dxw{x)pn{x)pm{x) = 6„ 


(18) 


Such polynomials can be generated using the recurrence 
[35] 


by diagonalizing the N x N matrix M [37] 


M = 


/ ao y/Pi 0 ... ^ 

v(3i ai y/132 
0 y/]^ ct2 




( 21 ) 


In addition, denoting the n-th eigenvector of M as the 
Christoffel weights in eq. (17) are given by the square of 
its first element 


Wn = vl^. ( 22 ) 

If the inner product (18) is not normalized, one has to 
multiply the right-hand side of this equation with the 
norm ^^dxw{x). 


2. Discrete Hamiltonian representation 

Let us now discuss in more detail how to obtain a dis¬ 
crete Hamiltonian with Nt, bath sites from the N roots 
Xm and Christoffel weights Wn that appear in the Gaus¬ 
sian quadrature rule for the integral (15). We discuss 
two cases (a) w{x) = J{x) and (b) w{x) = 1. Case (a) is, 
to our knowledge, the only one considered in the litera¬ 
ture [ 20 , 22-28], whereas case (b) makes the most simple 
choice for the weight function. 

(a) The choice w{x) = J{x) and fzix) = ^ 3 ^ leads 
to polynomials that are orthogonal with respect to 
J{x), which we therefore call bath-spectral-density- 
orthogonal (BSDO). Combining (15) and (17) it is 
found 


w 

A(z)«^—^=A<^-(z), (23) 

which allows to identify the Christoffel weights 
computed via ( 22 ) with the weights \Vn\^ of the 
discrete bath degrees of freedom 

1K^1' = Wn. (24) 


Pn+l(x) = (X- an)Pn(x) - /3nPn-l(x), (19) 

Po(x) = 1, p-iix) = 0, n = 0, ...,fV - 1, 

where /3o = 0 and 

7 „ = / dxpn{x)w{x), (20a) 

J a 

1 

an = — / dx x Pn{x)w{x), n = 0,...,N—l (20b) 

7n Ja 

I3n=lnhn-l, U = 1, N - 1. (20c) 

It is easy to see [36] that the roots of pn can be obtained 


(b) The choice w{x) = I and fz{x) = This is the 
case of Legendre polynomials and one obtains 


A(z) 


Nb 

yy 


WnJ{Xn) _ ^discr 




(25) 


n—1 


and the Christoffel weights Wn relate to the weights 
of the discrete bath via \ Vn\^ = WnJ{Xn). 

The next question is, which of these cases leads to a 
better approximation? Equations (23) and (25) derived 
from (17) do not hold exactly: in both cases (a) and (b) 
fz (x) contains a pole and hence it can not be exactly 


Z — X 
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represented by a polynomial of degree 2Nb — 1. Indeed, a 
pole is highly difficult to approximate with polynomials 
and it is quite irrelevant, whether one has an additional 
factor J{x) that multiplies this pole as in case (b), if 
this factor J{x) does not exhibit a severe non-regular 
behavior. This argument is confirmed by the numerical 
examples discussed in Section IV. 


where Unix) = yjJ{x)pnix). Therefore, they correspond 
to a weighted average over the total support of the spec¬ 
tral function J(a;). Note that due to orthogonality and 
normalization of p„(x), the transformation is unitary 

dxU*ix)Umix) = dxwix)pnix)pmix) = Snm and 
thereby retains the (anti-)commutation relation of aj,. 


3. Relationship to chain mappings 


4- Relationship to the Lanczos algorithm 


In this section, we show that the orthogonal polyno¬ 
mial method with the weight function chosen as wix) = 
Jix) (case (a) above), is equivalent to the chain mapping 
proposed in Refs. 20, 26, and 38, and recently modified 
in Ref. 28 to tackle temperature environments in an al¬ 
ternative way. It is also equivalent to the chain mapping 
derived in the Appendix of Ref. 24. The chain repre¬ 
sentation of the discrete star Hamiltonian obtained by 
considering wix) = Jix), can be written as 

= Hsys + Rtot id^eo + eld) (26) 

Ni-l Nb-2 

+ 'y ' Oinehen J y ) y/ fdn+lieh^ien -|- enCn+l), 

n—0 n—Q 


where |VtotP = f^^xJix) was defined in (12) and q;„ 
and Pn were defined in the recurrence relation (19). In 
the limit Nt ^ oo, becomes unitarily equivalent 

to the continuous H in (la), and thus provides an exact 
representation of H. 

For finite Nb, the unitary transformation that takes 
(26) back to its star representation (3), is equivalent to 
a diagonalization of the matrix (21) formed by the re¬ 
currence coefficients. As described above, such a trans¬ 
formation leads to the same weights and nodes as the 
ones obtained with the Gauss-Christoffel (BSDO quadra¬ 
ture). In other words, computing the system dynamics 
with a chain Hamiltonian (26) is equivalent to comput¬ 
ing the system dynamics with a star Hamiltonian (3) 
where nodes Xn and weights are computed with the 
BSDO quadrature. Regarding the important application 
of DMRG calculations: in contrast to what had been 
commonly believed, it was only recently shown that the 
star representation can be much less entangled than the 
chain representation [39]. 

Within the direct discretization strategy, the creation 
operators c^ of the discrete Hamiltonian in the star geom¬ 
etry (3) were obtained as an average over the continuous 
bath degrees of freedom aj, in a small interval as de¬ 
fined in (11). Within the orthogonal polynomial strategy 
described in the current section, the discrete operators in 
the chain Hamiltonian (26) are related to the continuous 
operators via 


/■b 

ei= dxUnix)al, 


(27) 


The measure wix) = Jix) is commonly known as Stilt- 
jes measure, and the three-term recursion (19) of the as¬ 
sociated BSDO polynomials is equivalent to the Lanczos 
algorithm for the continuous bath Hamiltonian i7bath in 
(la) [35] (see appendix G). The environment discretiza¬ 
tion then is a consequence of truncating the infinite re¬ 
currence relation (and therefore the matrix (21)) at a 
finite N = Nb- The implementation of the algorithm on 
a computer is though impossible, as there is no direct 
matrix representation for the continuous Hbath- 

By contrast, the Lanczos algorithm is a standard pro¬ 
cedure to tridiagonalize a given discrete bath Hamilto¬ 
nian Hbati[ (3) I obtain its unitarily equivalent 

chain representation. In order to so, one has to come 
up with a discrete Hamiltonian in the first place, which 
then has to be constructed using a direct discretization 
strategy. 


III. TIME EVOLUTION 


Let us now study the time evolution of the hybridiza¬ 
tion function, which describes the time evolution of the 
bath, and the time evolution of the Green’s function of 
the system, from which we can construct the time evo¬ 
lution of all system observables. The Green’s function is 
given by 


Git) = -f('(/'o|e 



\^o)=d^\Eo) 


(28) 


where the initial state is the excitation of the system 
iLsys through occupation with a particle, and the spectral 
density of the system is 

Aix) = ^ \iMEn)\Hix - iEn - Ao)), (29) 


where the sum is over all eigenstates \En) and eigenen- 
ergies of the full Hamiltonian (la). For a quadratic 
(single-particle) Hamiltonian, without loss of generality, 
one can consider Eq = 0 and \Eo) = Jvac) and therefore 
only has to study the time-evolution of a single particle 
that is initially in the system and starts interacting with 
the bath at non-zero times. 
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TABLE 1. Lanczos algorithm and orthogonal-polynomial strategy for real-time evolution. 


Lanczos algorithm 


Quadratic Hsya 

Non-quadratic Hays 

For continuous Hbath 

(lb)) 

(eq. 

"Hdiscr (40)) is obtained formally 

(App. Cl), and numerically (Sec. 

IIB 4) 

Same as for quadratic Hays 

For continuous H 
(la)) 

(eq. 

Hn (eq. (38)) is obtained formally 
(App. C 1 and Sec. Ill B 1 for first 
steps of algorithm). 

Not possible 

Is Lanczos for H equal to 
Lanczos for Hbath? 

Sec. IIIB2: Yes, Hn = for or¬ 

thogonal polynomial strategy 

Sec. Ill C: No, in general Hn ^ 


Analogously to (28), we define the time evolution of 
the hybridization function as 

/ OO 

dx J(a;)e““‘. (30) 

-OO 

For a discrete Hamiltonian ijdiscr^ obtains 

/ OO 

dxA'^‘"“(x)e-“*, (31) 

-OO 

Nb 

Adiscr(^) ^ |Kpe-“"‘. (32) 

n—1 

In the following, it is shown that the orthogonal poly¬ 
nomial strategy yields the best description of the short- 
and intermediate-time evolution of the continuous Hamil¬ 
tonian, if the latter is quadratic. It will then become clear 
why none of the discretization strategies can be consid¬ 
ered the best or the optimal one if the Hamiltonian is 
non-quadratic (has higher order interactions). In partic¬ 
ular: 

• Sec. HI A shows that the best approximation of 
(30) is obtained using the orthogonal polynomial 
strategy as described in Sec. HB. 

• Sec. HIBl shows that the Lanzos algorithm for 
the full H generates a matrix 'Htv, which gives 
the nodes and the weights that approximates the 
Green’s function (28) with a polynomial quadra¬ 
ture rule. 

• Sec. HIB2 shows that if idgys is quadratic, 'Hn = 
■^discr^ where T^^iscr jg obtained by Lanczos tridi- 
agonalization of idbath- Also, as it was shown in 
Sec. H B 4, a Lanczos tridiagonalization of Lfbath is 
equivalent to a bath discretization using the orthog¬ 
onal polynomial strategy of Sec.HB. Hence, the or¬ 
thogonal polynomial strategy leads to a quadrature 
rule also for the Green’s function (28). 

• Sec. HIC shows that if iLgys is non-quadratic, 
then 'Hn ^ and nothing can be concluded 

about the optimality of any particular discretiza¬ 
tion method. 

An overview of these steps is provided in Table (I). 


A. Time evolution of the bath 

In Sec. HB, we learned that polynomial quadrature 
rules provide us with the highest degree of exactness for 
computing the integral (17). In the following, we will 
see that this also helps us to understand in which cases 
(32) provides a good approximation of the Fourier type 
integral such as (30), and how to choose the parameters 
of the bath in order to obtain the best approximation. To 
this end, let us define the error term i?Vt(t) and write 

/ OO 

dx J(x)e-“* = ^ + i?v,(t). 

(33) 

We see that if we set w{x) = J(x) to construct orthog¬ 
onal polynomials via (19) and choose x„ to be the roots 
of the degree polynomial and = Wn to be the 
Ghristoffel weights (22), then (33) has the form of a Gaus¬ 
sian quadrature rule as in (17) with f{x,z = t) = 

That is, only if we choose \Vn\^ and x„ according to 
the orthogonal polynomial strategy with w{x) = J(x), 
our discrete Hamiltonian corresponds to evaluating the 
Fourier transform (33) to degree of exactness 2N{, — 1. 
Otherwise, the degree of exactness will be lower. What 
does this mean in practice? 

For a fixed time t, let us expand the part of the 
integrand J(x)e““* = iu(x)e““* in (33) that cannot be 
absorbed in a weight function in orthogonal polynomials 
q„(x), which are orthogonal with respect to v(x) (v(x) > 
0 is an arbitrary weight function), according to 

N OO 

^ ^ ^ Cn^n(^) H” ^ ^ ^nQn(^) ^ 

n—0 n—N+1 

Cn = f dxz;(x)e““‘g„(x). (34) 

J a 

Let us furthermore assume the family of polynomials 
qn{x) to be chosen optimally for the fixed time t. The 
optimal choice generates the most quickly converging 
sequence c„ —>■ 0 and by that minimizes the remain¬ 
der Tn = Y.n=N+lCn Pnix) at each order of N. Of 
course, we don’t know which polynomials these are, but 
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this is not relevant. The only property we need is that 
the coefficients become zero for values high values of n: 
c„ ~ 0 for n > N'{t), where N'{t) = i(6 — a)t (this is 
shown in Appendix D). 

The important observation to make is that choosing 
Xn and \Vn\^ = Wn according to the orthogonal poly¬ 
nomial strategy of Sec. IIB, corresponds to integrating 
the first term with N = 2Nh — 1 in (34) exactly. Any 
other choice, will lead to an exact integration of the term 
only at a lower order, or will not integrate it exactly at 
any order. Combining this observation with the fact that 
c„ ~ 0 for n > ^{b—a)t, we conclude that the orthogonal 
polynomial strategy reproduces basically the exact time 
evolution of the hybridization function for t < tmax, with 


Assume iJsys = Socf'd quadratic. Let us take as initial 
Lanczos vector the state |/o) = |d) = d'^|vac) = |^/>o)- 
Denoting the single-particle states of the bath as la^;) = 
aj,|vac), we have following (Cl) 

ao = (/o|-ff|/o) = £0) 

\ro) = H\fo) - aolfo) = f dxV{x)\ax), 

J a 

(^-oko) = dx\V{x)\'^ = |VtotP = Pf, 

J a 

l/i) = ^ / dxV{x)\aP). (37) 

'^tot Ja 


t 


max 


2 ?^ 
b — a 


Continuing the algorithm up to order N produces a trun- 
cated representation of H, which is a, N x N matrix. 


This result is confirmed in the numerical experiments in 
Sec. IV. We have therefore shown that the best approxi¬ 
mation of (30) is given by a orthogonal polynomial strat¬ 
egy as described in Sec. IIB. 


B. Time evolution of the system 

The Green’s function of the system as defined in (28) 
can be rewritten as follows 

/ OO ^ 

dxA(x)e-“* = ^ |(V'o|A„_i)pe-*^--‘ 

n -1 

N 

= Y,\{Wn)\^e-^^-*+RN{t), (36) 

n—1 

where |i?„) are eigenstates and eigenenergies of the 
exact, continuous Hamiltonian (la), and RN(t) is a re¬ 
mainder. The problem is therefore again to choose the 
states \Xn) and the nodes A„, such as to make (36) a 
quadrature rule, which we just showed (Sec. Ill A) to 
yield the best approximation of Fourier type integrals. 


1. Lanczos for quadratic Hamiltonian 



/ £o 

Vtot 

0 


\ 



Vtot 

5i 


0 


Rn = 

0 


Ct2 



■ (38) 


0 

0 


0-3 




1 : 




J 


As discussed in 

Appendix C 

, there is a 

set of orthogo- 


nal polynomials qn(x) that are orthogonal with respect 
to w{x) — A{x) (A(x) is the spectral density of the 
full Hamiltonian H) associated with the preceding Lanc¬ 
zos algorithm. Therefore, diagonalization of (38) yields 
roots Xn and Christoffel weights Wn = |(/o|V„)p = 
|(^/>o|Ar„)p. Hence, the Lanczos algorithm evaluated 
for the continuous quadratic Hamiltonian H with initial 
state l/o) = \ipo) generates the nodes and weights that 
make the approximation (36) a quadrature rule. Note 
that Xn ^ En, since En-i are true eigenvalues of H, 
and Xn are the eigenvalues of the truncated tri-diagonal 
representation Rjq of H. 

But how does this relate to the parametrization for a 
discrete Hamiltonian that we obtain from the or¬ 

thogonal polynomial strategy (19) for the weight function 
w(x) = J(x)7 


For quadratic Hamiltonians we will show in the follow¬ 
ing, that the orthogonal polynomial strategy (19) gener¬ 
ates a quadrature rule for (36), and A„ and \Xn) become 
respectively the eigenenergies and eigenstates of the dis¬ 
crete Hamiltonian i/^iscr does not use the 

orthogonal polynomial strategy, or the Hamitonian is not 
quadratic, one never generates a quadrature rule in (36). 

To this end, let us compute the first steps of the stan¬ 
dard Lanczos tridiagonalization algorithm recapitulated 
in Appendix C. Here, we do it for the full continuous 
quadratic Hamiltonian (la), and not for the bath and 
coupling part of the discrete Hamiltonian (3), as usually 
done in the context of chain mappings. 


2. Equivalence with orthogonal polynomial strategy 

The discrete quadratic Hamiltonian iJ^iscr^ which has 
dimension (Nb -I- 1 ) x (Ni, + 1 ), generates the following ap¬ 
proximation to the time evolution of the Green’s function 
of the continuous system 

/ OO 

da;A(a;)e-“‘ 

-OO 

iVb + l 

= E * + <---(t), (39) 

n—1 





where are eigenstates and eivenvalues of 

^discr ^discr |^g represented in the chain ge¬ 

ometry (26) as 


/^discr _ 


/ 


£0 

f4ot 

0 

0 


f4ot 0_ ... \ 

0^0 V/^i 0 

\/Wi 0^1 's/S ' ■ 

0 y/h 0^2 


V : ■•• ■■• / 


(40) 


As l^o) = M), this representation of iJ^iscr directly yields 
the weights and energies in (39). 

In the following, we will show that the matrix (40) 
equals the matrix (38) that generates the quadrature 
rule, only if we compute the parameters of the discrete 
Hamiltonian using the orthogonal polynomial strategy 
(19) with w{x) = J{x). Only then, also (39) is a quadra¬ 
ture rule. 

To this end, let us further evaluate the Lanczos algo¬ 
rithm for the continuous H. Using the results of (37), 
we can represent the terms in (la) as 


Hsys — £o|/o)(/o|, 

-ffcoupl = Kot(|/o)(/l| + h-C.), 

77 bath = / dxx\aa:){ax\. 

J a 

As the Lanczos basis is orthogonal, we see that in subse¬ 
quent Lanczos steps, only iLbath can contribute: i7sys and 
77coupi only have contributions in the subspace spanned 
by l/o) and |/i). We therefore have to evaluate a single 
next Lanczos step using the full H, and from then on 
can iterate using only Llbath- Now note that the Lanc¬ 
zos vector l/i) in (37), which is the starting vector for 
subsequent Lanczos steps, equals the state |eo) in (C4), 
which is the initial state for a tridiagonalization of the 
bath. We already know the latter to be equivalent to the 
orthogonal polynomial strategy. The Lanczos recursion 
for the full H therefore generates the coefficients of the 
orthogonal polynomial strategy. Let us check this for the 
next step. 


5 l = (/l|i 7 |/l) = (/l|i 7 bath|/l). 

In) = H\fi) - 3i|/i) - Utotl/o) 

= 77 bath|/l) — ai|/i). 


Evidently, 5i = ao and |ri) = |ro) as |/i) = |eo) such 
that this equals the parameters of (C5) and (19). Hence 
the matrices (40) and (38) are equivalent, and the time 
evolution computed with the discrete Hamiltonian is a 
quadrature rule. 

For any other choice of which is not 

parametrized using (19), we would not obtain an equiva¬ 
lent representation to (38), and therefore, (39) would not 
be a quadrature rule. 


The estimate (35) for the maximal time tmax yields, 
as the quadrature rule now uses a polynomial of degree 
A’b + 1) 


t 


max 


^ 2Nb + l 
b — a 


(41) 


C. Impossibility of optimal choice for 
non-quadratic Hamiltonians 

If Hsys is not quadratic, but has higher order in¬ 
teraction terms, we cannot obtain a representation of 
^discr terms of single-particle states, and hence as 
a {Nb -I- 1) X {Nb + 1) matrix. Rather, any repre¬ 
sentation of then has an exponential dimension, 

e.g. 2'^*’+^ X 2^*’+^ for spinless fermions, and dimension 
jjNi+i ^ jjNi +1 tQj. bosons with a local basis truncated 
at a dimension D. The summation over the discrete time 
evolution of (39) then involves an exponential number of 
terms. By dimensionality, this summation can never cor¬ 
respond to a quadrature rule with Nb parameters, which 
gives rise to Nb roots. The time evolution of the bath 
hybridization function, which always is a single-particle 
evolution, is not affected by this argument and is still 
best described using the parameters provided by the or¬ 
thogonal polynomial strategy. 

In summary, for non-quadratic Hamiltonians, even if 
we have a good approximation of the bath hybridization 
function up to tmax, the dynamics of the system, given 
by the Green function (39) will no longer be exact up to 
this time. 


IV. NUMERICAL EXAMPLES 
A. Spin-boson model 

Let us consider the Hamiltonian of an OQS with Hsys 
coupled to a continuous bosonic reservoir 


H = Hs. 


dkg{k){b{k)a^+a b{k)^) 


dk U!{k)b{ky b{k), 


(42) 


where g{k) are the coupling strengths, and b{k) (5(fc)^) 
are harmonic oscillator operators with commutation re¬ 
lations [b{k),b{k')^ = S{k — k'). Here, the index k labels 
the modes, which have a maximum momentum fcmax- In 
the frequency representation, and provided that the envi¬ 
ronment is initially in a Gaussian state, this Hamiltonian 
can be rewritten as 


H = Hs. 


dujg{(jj) + b(uj)^a ) 


du) b{Lo)^b{w), 


(43) 
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where o-’m ax is determined by fcmax, and we have defined 
g{uj) = J{uj), where J{uj) = 5 ^(w)pdos(w) is the spec¬ 
tral density of the environment, and Pdos(w) is the envi¬ 
ronment density of states. Hence, the Hamiltonian (43) 
acquires the form (la), obviously once interpreting the 
continuous variable x as uj, and d = a~. We also note 
that the above Hamiltonian corresponds to a simplified 
version of the spin-boson model, as it assumes a rotat¬ 
ing wave approximation to discard fast rotating terms of 
the form b^{k)(7'^, and b{k)a~. Such an approximation, 
which is particularly valid in quantum optics, leads to a 
Hamiltonian that conserves the number of particles. This 
simplifies considerably the numerical treatment, particu¬ 
larly at zero temperature. 

In order to characterize the environment, let us con¬ 
sider a spectral density of the Caldeira and Leggett type 
[40, 41], 

J(w) (44) 

which constitute a very general description that allows 
to describe many different types of reservoirs, depend¬ 
ing on the choice of the parameter s. The exponential 
factor in this model provides a smooth regularization for 
the spectral density, being modulated by the frequency 
iOc- Environments with 0 < s < 1 are considered as sub- 
ohmic^ while those corresponding to s = 1 and s > 1 are 
known as ohmic and super-ohmic respectively. The con¬ 
stant a describes the coupling strength of the system and 
the environment. In the following, we will focus on a sub- 
ohmic spectral density with s = 1/2. Sub-ohmic spectral 
densities describe the frequency dependence of photonic 
bands in photonic band gap materials [27, 42, 43], as 
well as the dominant noise sources in solid state devices 
at low temperatures such as superconducting qubits [44] , 
nanomechanical oscillators [45], and quantum dots [46]. 

Considering zero temperature, the OQS dynamics can 
be easily solved by exact diagonalization (ED), since 
there is only one excitation involved in the problem (it is 
a single-particle problem with a quadratic Hamiltonian). 
In this context. Fig. (1) shows results for the population 

P{t) = {a+{f)a-{t)) (45) 

£:(t) = |P(t)-P‘i-“(t)|, (46) 

where P(t) is computed with the continuous environ¬ 
ment, and is the population computed with the 

discretized environment. Eft) is the error made by us¬ 
ing the discretized environment. We compare results ob¬ 
tained using the linear discretization as an example for 
a direct discretization strategy with the orthogonal poly¬ 
nomial strategy that uses (19) with the weight function 
w{x) = J(x) generating BSDO polynomials. Clearly, the 
BSDO strategy leads to an error that is at least two or¬ 
ders of magnitude smaller than the one of the linear dis¬ 
cretization with the same number of modes until reach¬ 
ing a time ttmax, when the discretized system fails to 
accurately describe the continuous system. Physically, 
such a failure can be interpreted as a revival of the sys¬ 
tem dynamics, which occurs when the emitted excitation 



t 


FIG. 1. Time evolution of the population of the upper level 
for Ni, — 65 (upper panel) and error £ according to (46) 
in logarithmic scale (lower panel). In both cases, different 
discretization schemes are considered. Dot-dashed green and 
dashed orange curves correspond respectively to polynomial 
and linear methods. The linear black curve in the upper panel 
corresponds to the exact solution. The curves and error of the 
mean and the equal weight method of Sect. (IIA 1) are not 
shown, but have a similar behaviour as the ones of the linear 
method. The red line below shows the time tmax at which 
the error of the polynomial method increases two orders of 
magnitude,which coincides with the exact formula (41) (see 
also Fig. (3)). We have considered Ua = 0.5, a = 1, s = 0.5, 
ojc = 10, and a maximum frequency in the spectrum tUmax = 
50. 


hits the chain extreme and bounces back into the sys¬ 
tem. We note that the results obtained with the heuris¬ 
tic approaches described in Sec. (HA 1) (not shown), are 
found to achieve a similar level of accuracy as the linear 
discretization strategy. 

Fig. (2) compares two orthogonal-polynomial based 
strategies: one generated with (19) using the weight func¬ 
tion w{x) = J{x) (BSDO quadrature) and one using 
w{x) = 1 (Legendre quadrature). The figure confirms the 
statement made after eq. (25) that both strategies yield 
basically the same accuracy if the bath spectral density 
does not show a severe non-regular behavior. Also, as 
shown in Fig. (3), tmax is linearly related to the number 
of node points considered in the quadrature rule. This 
follows from equation (41). 

We note that also in the finite temperature case, stud¬ 
ied within the second order weak coupling master equa¬ 
tion, allows us to recover the result that the BSDO strat¬ 
egy is optimal up to the time tmax (see Appendix (E)). 
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Exact 

A 

Linear 

& 

0 

Mean 

■ 

BSDO poly 



FIG. 2. Evolution of the upper level considering the quadra¬ 
ture method with different polynomial classes for Nb — 65 
nodes. Blue diamonds, and green squares correspond, respec¬ 
tively, to the Gaussian quadrature rule (with Legendre poly¬ 
nomials), and to the Gauss-Christoffel quadrature with BSDO 
polynomials (i.e. polynomials obeying the relation (18) with 
w{x) = J{x)). 



FIG. 3. Maximum time at which the error between the evolu¬ 
tion with discretization with N nodes, and the exact (contin¬ 
uous) one is below a certain threshold chosen as 0.004. The 
maximum frequency in the spectrum is tUmax = 100. Blue 
diamonds, and green circles correspond, respectively, to the 
Gaussian quadrature rule (with Legendre polynomials), and 
to the Gauss-Ghristoffel quadrature (with BSDO polynomi¬ 
als). System parameters are the same as in Fig. (1), except 
for the fact that we are now considering s = 1.5. The figure 
shows approximately the same slope as the one predicted by 
eq. (41). 


B. Single-impurity Anderson model 


The single-impurity Anderson model (SIAM) has the 
form of Hamiltonian (la), with the impurity and bath 
operators being spin-dependent fermionic creation and 


FIG. 4. The generic bath spectral density (48) and its 
discretized versions. To plot the discrete spectral function 
jdiscr(^), replace the delta function by a rescaled indica¬ 
tor function 5{x — x„) —^ x(® “ x„)/Aa;„, where Axn is the 
width of In = [{Xn + Xn-i)/‘2, (xn+i + Xn)l‘2\- This rescaling 
accounts for the fact that for comparisons with the contin¬ 
uous spectral density, the discrete spectral function should 
be interpreted as a probability density defined on the energy 
interval (a, b) that associates a weight (an excitation proba¬ 
bility) to an energy interval, and not as a probability mass 
function that associates a weight to a value of Xn- 


annihilation operators, 

Hsys = U{d)^d'f — ~ 2 ^’ 

.^bath — ^ ^ / dxXCi^^Q^xa-! 

a 

pb 

.^coupl - E dx V{x)(dl 

^X(T -bh.c..) 

^ a 


(47) 


In a grand-canonical picture this corresponds to the half- 
filled case obtained for chemical potential fi = —17/2. 
The physics of this case shows generic features. Clearly, 
for [/ yf 0, Hsys describes a non-quadratic interaction. 

The generic case of interest for the physics of strongly- 
correlated electron systems is best captured by a bath 
spectral density of the form 

= E^o 6 {- 4 ,o, 4 } e ^ for X G [-5,5] (48) 

outside of the interval [—5,5] we set J{x) — 0. This 
bath spectral density is a superposition of three Gaus¬ 
sian peaks that produces “gapped” regions where J{x) is 
practically zero. Figure 4 shows the continuous and the 
discretized version of this J{x). The peak at zero fre¬ 
quency corresponds to low-energy excitations in the bath, 
as they are present in a metal. The two other peaks cor¬ 
respond to high-energy excitations that become relevant 
when the interaction U generates low (single occupation) 
and high (double or zero occupation) energy states. In a 
Mott insulator, there is no low energy physics any more 
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FIG. 5. Time evolution of the SIAM (47) for [7 = 0. Upper 
panel: Time evolution for Nb = 15. Lower panel: Error for 
Nb = 31. The maximal time (red vertical line) until which the 
BSDO polynomial discretization yields the exact description 
can be computed using (41), and yields for a = —5, b = 5 and 
Nb = 31 the value tmax = 12.6. 



t 


FIG. 6. Time evolution of the SIAM (47) for [7 = 4. Upper 
panel: Time evolution for Nb = 15. Lower panel: Error for 
Nb = 31. 


and the interaction created a gap in the excitation spec¬ 
trum. The most exciting physics happens in the inter¬ 
mediate regime where the quantum Mott-Insulator phase 
transition occurs. 

Let us first study the non-interacting case [7 = 0, 
which only involves a quadratic Hamiltonian. In this 
case, we confirm the results of the previous section. Fig¬ 
ure 5 shows the time evolution of the overlap of the initial 
state (the Green’s function iG{t) = 
defined in (28)), that consists in placing a spin-up elec¬ 
tron on the impurity {ipit = 0)) = (i||i7o), with its time 
evolution. Evidently, the linear discretization yields the 


worst results, and the Gauss-Christoffel (BSDO) strategy 
yields a numerically exact result up to time 6. 

Let us now turn to the interacting case where U is non¬ 
zero and the Hamiltonian is no longer quadratic. Figure 
6 confirms the result of Sec. HI that BSDO polynomi¬ 
als do no longer give optimal results as they no longer 
generate a Gaussian quadrature rule. Now the heuristic 
mean method produces the best results, leading to errors 
that are at least a factor 2 smaller than the BSDO strat¬ 
egy. The mean method directly uses the fact that one 
can ignore gapped regions in the bath spectral density. 
This is important in the computation of strongly corre¬ 
lated materials. In both cases described in Figs. 5 and 6, 
the equal weight method of Sec. HA 1 performs qualita¬ 
tively similar to the mean method, and therefore it has 
not been shown for the shake of clarity in the figure. 


V. CONCLUSIONS 

In this paper we have analyzed a OQS coupled to 
a bosonic environment characterized by a Galdeira and 
Leggett type of spectral density, and a quantum impurity 
model consisting on an impurity coupled to a fermionic 
bath. We considered direct discretization strategies and 
orthogonal polynomial quadrature based strategies. We 
have shown that when using orthogonal polynomials, the 
choice of the polynomial class does not affect consider¬ 
ably the error in the resulting system dynamics. In addi¬ 
tion, we have shown that the Gauss-Christoffel quadra¬ 
ture rule (which is based on the choice of a particular fam¬ 
ily of polynomials here denoted as BSDO), correspond to 
the chain mapping approach proposed by Refs. 26 and 47. 
Such chain mapping leads effectively to a discrete chain 
representation, which when transformed back to a diago¬ 
nal form, leads to environment eigenvalues that precisely 
correspond to the nodes of the Gauss-Christoffel (BSDO) 
quadrature rule. 

Finally, we have shown that in a non-interacting sys¬ 
tem (i.e. with quadratic Hamiltonian), the polynomial 
quadrature method is exact at short times. Neverthe¬ 
less, for non-quadratic Hamiltonians (like an an impurity 
with non-zero interaction term) this is no longer the case. 
This means that the notion of optimality that is associ¬ 
ated with an optimal representation of the continuous in¬ 
tegral of J {x) by a finite number of points breaks down if 
we consider non-quadratic Hamiltonians. In other words, 
the non-linear problem that is encoded in such non¬ 
quadratic Hamiltonian obviously will no longer be well 
described by just considering a polynomial quadrature 
rule on the integral. It is noted that, although we have 
presented a scheme (the mean method) that performs 
better than Gauss Christoffel quadrature in this case, we 
showed that a general statement cannot be made. 

Note added in proof. Dynamical error bounds on ex¬ 
pectation values of system observables for a Hamiltonian 
discretised using orthogonal polynomials have recently 
been derived in Ref. [48]. 
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extend the promising DMRG calculations for situations 
with a non-entangled initial state [52? ] to the general 
case of entangled initial states. But then one also has to 
discretize the “first” bath, which incorporates the spec¬ 
tral information of H and which is equivalent to the bath 
that is the subject of the present paper. For the first 
bath, one again faces the problem that a cost function 
cannot be meaningfully defined. 

Appendix B: System Green’s function 


Appendix A: Numerical optimization 

Numerical optimization can be formulated in a 
straightforward way for the hybrdization function A(z), 
defined in (5), evaluated on a grid of imaginary frequen¬ 
cies z = iwk, 

x2 = ^|A(*a;fc)-Ad—(zccfcjp, (AI) 

k 

using standard numerical minimization techniques [31, 
49, 50]. On the real axis, the equivalent cost function can 
be formally defined as |A(w/c-|-jO''") —A'^’®“''(tLH;-|- 

but is of no use as the difference of a continu¬ 
ous function and a singular function is always infinite. 
Therefore, we cannot use numerical optimization to dis¬ 
cretize the spectral representation of the continuous bath, 
i.e. the hybridization function evaluated on the real axis 
via J{x) = —ilmA(a; -I- z0+). 

If one carries out the optimization on the imaginary 
axis via (Al), one obtains a set of parameters {x„, U„} for 
the discrete bath and an associated hybridization func¬ 
tion A‘i'"'’^(z), which gives a quantitatively precise ap¬ 
proximation to A(z) only when evaluated on the imag¬ 
inary frequency axis. On the real-frequency axis, the 
approximation is very rough and can only be considered 
qualitatively correct. This follows already from the fact 
that only relatively small numbers of bath sites Ab < 15 
can be stably optimized. Still the approach is valid if one 
is satisfied with the much lower precision on the real axis 
and does not strive to describe real-time evolution as in 
this paper. The preceding statements are e.g. discussed, 
among several other results, in Ref. 51, where the goal 
was not to describe real-time evolution but “thermody¬ 
namic” properties. 

We note that one can define a meaningful cost function 
on the real axis, if one allows for non-hermitian Hamil¬ 
tonians with complex bath energies, or an equivalent de¬ 
scription in terms of Lindbladt operators [32] . 

We further note that one can also construct an optimal 
discrete representation of the “second bath” that appears 
within non-equilibrium DMFT [33] . But this only suffices 
to describe situations in which the system and bath are 
initially not entangled [33] . As non-equilibrium DMFT is 
a promising approach to describe the non-equilibrium dy¬ 
namics of strongly correlated materials, it is desirable to 


The retarded system Green’s function is defined in 
terms of the general retarded Green’s function (system 
and bath) 


G(x) 


1 

X + iO — {H — Eq) 


(Bl) 


by taking expectation values with respect to the system 
states [19], e.g. l^/>o) = d^Eo) 


Gsys{x) = {tjjo\G{x)\tlJo). (B2) 

For the system Hamiltonian iLgys = cod^d it reads [19] it 
can be evaluated as 


^ + iO-eo+A(x)- 
where A(a;) is defined in (5). 


(B3) 


Appendix C: Lanczos algorithm 

1. General Lanczos algorithm and relation to 
orthogonal polynomials 

The Lanczos algorithm constructs a three-diagonal ma¬ 
trix representation of any Hermitian operator H by rep¬ 
resenting it in its Gram-Schmidt orthogonalized Krylov 
basis {1/n)}: Given a start vector j/g) that has non-zero 
overlap with all eigen-states of H, one orthogonalizes the 
vector \fn) with respect to all previous vectors ]/„') with 
n' < n. This results in 

an = {fn\H\fn), 

\rn) = H\fn) - an\fn) “ \/Kl/n-l) 

Pn+1 = \{rn\rn)\ , Po = 0 

\fn+i) =— 7 ^=\rn), for n = 0, ...,fV-1. (Gl) 
v Pn +1 

One can show that the Lanczos algorithm implicitly 
constructs a family of polynomials qn{x) that are orthog¬ 
onal with respect to an inner product weighted with the 
spectral density A{x) of the operator H [35, 53] 

dim(ff) 

w{x)= Y. \{En\h)\''5{E - E^) = A{x). 

n—1 
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The proof is as follows. Let us define the polynomial 
qn{x) of degree n via 

\fn)=qn{H)\h), (C2) 

and then show that they are orthogonal with respect to 
A{x). We note that (C2) can always be fulfilled as |/„) 
is constructed by applying H n times to the initial state 
l/o). Furthermore, 

pb 

/ dxA{x)qk{x)qi{x) = '^{fo\En)qkiEn)qi{En){En\fo) 

n=l 

= {fo\qkiH)qi{H)\fo) = ifklfi) = Ski, 
which completes the proof. 



2. Chain mapping 


In the following, we show how to use the Lanczos al¬ 
gorithm to tridiagonlize the star Hamiltonians E[ in (la) 
and in (3) . This amounts to using the general algo¬ 

rithm (Cl) for the bath Hamiltonians Ffbath and Hbath^, 
respectively. The bath Hamiltonians are quadratic and 
therefore simple to treat. They have the spectral den¬ 
sities J{x) and J‘^'^^''{x) as defined in (2) and (4), re¬ 
spectively. Already from this we can conclude from the 
argument of Sec. C 1, that the Lanczos algorithm applied 
for the continuous iLbath, yields the same set of orthog¬ 
onal polynomials as the recurrence (19), and is therefore 
equivalent to it. 

In practice, the algorithm is usually used to obtain rep¬ 
resentations of the discrete bath and coupling Hamiltoni¬ 
ans Llbat™ LIcou^i- We will lay out the procedure for 
the discrete case, and note differences to the continuous 
case where necessary. 

Let us denote the (single-particle) bath orbital states of 
the discrete star representation (3) as |c„). These are as¬ 
sociated with the operators via |c„) = cjjjvac). Anal¬ 
ogously, define the bath orbitals of the chain representa¬ 
tion (26) as |e„) where |e„) = ejj|vac). The first orbital 
of the chain representation then is 


JVb 


|eo) = ^^Vn\Cn), 
vt-^ 


(C3) 


n—1 


n—1 

in the discrete case, and 

rb 


|HotP= [ dxJ{x), 

r,-i -'a 


|eo) = [ dxV{x)\a:r), la^;) = 4|vac), (C4) 

Ltot Ja 


in the continuous case, in agreement with (27). In both 
cases, it is a superposition of all states in the star. The 
coupling Hamiltonians ™ (3) can then be written 

-^cou^i = ^ot(|d)(eo| -I-h.c.), where \d) is associated 


FIG. 7. Blue curves represent the master equation solution 
for the atomic population {t)a{t)), with quasi-continuous 
spectrum (plain solid curve) and Gauss-Christoffel (BSDO) 
quadrature with N = 100 nodes (curve with triangles). Or¬ 
ange curves (see also inset) represent the evolution of r(t) = 
fg aT('r)e*“®’^ for the same two cases. The spectral density, 
as well as all system parameters are the same as in Fig. (1), 
except for the coupling that now is considered to be weak, 
a = 0 . 01 . 


with the system operator rfl. The same equation holds 
in the continuous case. 

One then uses the Lanczos algorithm to construct a 
three-diagonal representation of i7bath^ 


a„ = (e„|i7btth |en), (C5) 

kn) = -ffb^th |en) - a„|e„) - V^|e„-i) 

/3„+i = |(r„|r„)|, /3o = 0 

1 


|en-i-i) — 




), for n = 0 ,..., A^b — 1. 


n+1 


or analogously, for the continuous case. The parameters 
a„ and /3„ in the recursion are the parameters of the 
Hamiltonian (26), and with that the map is complete. 

In practice we note that we cannot find a direct matrix 
representation of the continuous Hamiltonian (la) that 
we could use on a computer to compute (C5). In the 
discrete case, on the other hand, the preceding equations 
are easily solved by generating a matrix representation by 
multiplying from the left with {cn> \ and inserting identi¬ 
ties |cn') (cn' | such that the initial vector can be writ¬ 
ten as ((c„|eo))()Ci = (14)^ii the representation of 
^bith involved is (c„|i7bath|c„') = 

Due to the numerical instability of the Lanczos algo¬ 
rithm, the recurrences (C5) and (19) have to be computed 
with high-precision arithmetics when exceeding Ni, ~ 40 
or using the stabilized implementation of Ref. 35. 
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FIG. 8. Same as in Fig. (7), but considering finite tempera¬ 
ture (/3 = 1). 


Appendix D: Estimate the error in time evolution 

As Chebyshev polynomials are almost optimal they 
will result in a sequence c„, which is very close to the 
sequence produced by an optimal choice of polynomials, 
in the sense of the discussion of (34). 

For Chebyshev polynomials {v{x) = v{x') = ^(1 — 
x')~i and qn{x) = qnix') = arccos(ncos(a;')) with x' = 
— I, X = ^{b—a)x' + ^{b + a)), we can evaluate the 
coefficients in (34) explicitely, 

= (Dl) 

where J„(t') are Bessel functions of the first kind. For all 
practical purposes, ~ 0 if n > t'. More concretely, 

the asymptotic form for high values of n reads n ^ — 1, 


Jn{t') ^ (n+i)! (t)” shows that this decreases 

as a faculty. 

Appendix E: Open quantum system in the presence 
of a thermal environment 

Let us now consider d = = ax in (la), and a h- 

nite temperature in the environment. We study this 
case within a standard approximate scheme, namely a 
master equation (ME) up to second order in the system- 
environment coupling parameter g [ 1 ], 

^^df^ = -i[Hsys, Psit)] + J dra^it - T)[d\ ps{t)d{T - t)] 

+ [ dTa2{t-T)[d\T-t)psit),d] 

Jo 

+ I draiit - T)[d{T - t)ps{t),d'^] 

Jo 

+ [ dTal(t-T)[d,Ps{t)d{T-t)'^]+0{g^), (El) 
Jo 

with ai{t -t) = Y^kdlink + a 2 {t - t) = 

and d{t) = 

As it can be seen in Figs. 7, for zero temperature, and 8 
for finite temperature, the polynomial Gauss-Christoffel 
(BSDO) quadrature is still extremely accurate a short 
times. Nevertheless, just as in the zero temperature 
case, after a certain time tmax, the discretization pro¬ 
cedure starts to fail. Such a failure is originated from 
the fact that the polynomial quadrature rule starts to re¬ 
produce inaccurately the integrals T{t) = aT(T)e*‘^®'^, 

with arit) = ai{t) + cx^it), entering in the master equa¬ 
tion. 

Indeed, as seen in the inset of both figures, small devia¬ 
tions of this quantity due to an inaccurate discretization, 
produce large deviations in the dynamics with respect to 
the reference (corresponding to the solution with a quasi- 
continuous spectrum), and this deviation is particularly 
large at hnite temperatures. 


[1] H. Breuer and F. Petruccione, The theory of Quantum 
Open Systems (Oxford Univ. Press, 2002). 

[2] A. Rivas and S. F. Huelga, Open Qnantum Systems. An 
Introduction (Springer, Heidelberg) (2011). 

[3] P. W. Anderson, Phys. Rev. 124, 41 (1961). 

[4] W. Metzner and D. Vollhardt, Physical Review Letters 
62, 324 (1989). 

[5] A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992). 

[6] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozen- 
berg. Rev. Mod. Phys. 68 , 13 (1996). 

[7] G. Kotliar, S. Savrasov, K. Haule, V. Oudovenko, O. Par- 
collet, and G. Marianetti, Reviews of Modern Physics 78, 
865 (2006). 

[8] T. Maier, M. Jarrell, T. Prnschke, and M. Hettler, Rev. 


Mod. Phys. 77, 1027 (2005). 

[9] D. Zgid and G. K.-L. Chan, The Journal of Chemical 
Physics 134, 094115 (2010), arXiv: 1012.3609. 

[10] W. Domcke, Physics Reports 208, 97 (1991). 

[11] F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, 
Phys. Rev. B 79, 235336 (2009), arXiv;0903.2414. 

[12] R. Peters, Phys. Rev. B 84, 075139 (2011). 

[13] F. Guettge, F. B. Anders, U. Schollwoeck, E. Eidel- 
stein, and A. Schiller, Phys. Rev. B 87, 115115 (2012), 
arXiv: 1206.2186. 

[14] M. Ganahl, P. Thunstrom, F. Verstraete, K. Held, and 
H. G. Evertz, Phys. Rev. B 90, 045144 (2014). 

[15] F. A. Wolf, I. P. McCulloch, O. Parcollet, and 
U. Schollwock, Phys. Rev. B 90, 115124 (2014), 















15 


arXiv:1407.1622. 

[16] A. Weichselbaum, F. Verstraete, U. Schollwock, J. Cirac, 
and J. von Delft, Phys. Rev. B 80, 165117 (2009). 

[17] F. A. Wolf, J. A. Justiniano, I. P. McCulloch, and 
U. Schollwock, Phys. Rev. B 91, 115144 (2015), 
arXiv:1501.07216. 

[18] M. Zwolak, The Journal of Chemical Physics 129, 101101 
(2008). 

[19] A. C. Hewson, The Kondo Problem to Heavy Fermions 
(Cambridge University Press, 2003). 

[20] A. W. Chin, S. F. Huelga, and M. B. Plenio, Semiconduc¬ 
tors and Semimetals, edited by U. Wurfel, M. Thorwart, 
E. R. Weber, and C. Jagadish (Academic Press, 2011). 

[21] R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 
80, 395 (2008). 

[22] R. S. Burkey and C. D. Cantrell, J. Opt. Soc. Am. B 1, 
169 (1984). 

[23] A. K. Kazansky, Journal of Physics B: Atomic, Molecular 
and Optical Physics 30, 1401 (1997). 

[24] M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B 77 , 
075116 (2008). 

[25] N. Shenvi, J. R. Schmidt, S. T. Edwards, and J. C. Tully, 
Phys. Rev. A 78, 022502 (2008). 

[26] J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, 
Phys. Rev. Lett. 105, 050404 (2010). 

[27] J. Prior, I. de Vega, A. W. Chin, S. F. Huelga, and M. B. 
Plenio, Phys. Rev. A 87, 013428 (2013). 

[28] I. de Vega and M.-C. Banuls, (2015), arXiv:1504.07228. 

[29] F. A. Y. N. Schroder, A. W. Chin, and R. H. Friend, 
(2015), 1507.02202. 

[30] I. de Vega, Phys. Rev. A 90, 043806 (2014). 

[31] M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 
(1994). 

[32] A. Dorda, M. Nuss, W. von der Linden, and E. Arrigoni, 
Phys. Rev. B 89, 165105 (2014). 

[33] C. Gramsch, K. Balzer, M. Eckstein, and M. Kollar, 
Phys. Rev. B 88, 235106 (2013), arXiv; 1306.6315. 

[34] W. Gautschi, “A survey of gauss-christoffel quadrature 
formulae,” (Birkhauser Verlag, Basel, 1981) pp. 72-147. 

[35] W. Gautschi, Journal of Computational and Applied 
Mathematics 178, 215 (2005). 

[36] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and 


B. P. Flannery, Numerical Recipes 3rd Edition: The Art 
of Scientific Computing , 3rd ed. (Cambridge University 
Press, New York, NY, USA, 2007). 

[37] G. H. Golub and J. H. Welsch, Mathematics of Compu¬ 
tation 23, 221 (1969). 

[38] M. P. Woods, R. Groux, A. W. Chin, S. F. Huelga, 
and M. B. Plenio, Journal of Mathematical Physics 55, 
032101 (2014). 

[39] F. A. Wolf, 1. P. McCulloch, and U. Schollwock, Phys. 
Rev. B 90, 235131 (2014). 

[40] A. O. Caldeira and A. J. Leggett, Physica A 121, 587 
(1983). 

[41] U. Weiss, Quantum Dissipative Systems, edited by 
W. Scientihc (Series in modern condensed matter sys¬ 
tems, 2008). 

[42] M. Florescu and S. John, Phys. Rev. A 64, 033801 (2001). 

[43] I. de Vega, D. Alonso, and P. Gaspard, Phys. Rev. A 
71, 023812 (2005). 

[44] A. Shnirman, Y. Makhlin, and G. Schn, Physica Scripta 
2002, 147 (2002). 

[45] G. Seoanez, F. Guinea, and A. H. C. Neto, EPL (Euro¬ 
physics Letters) 78, 60002 (2007). 

[46] N.-H. Tong and M. Vojta, Phys. Rev. Lett. 97, 016802 
(2006). 

[47] A. W. Chin, J. Prior, S. F. Huelga, and M. B. Plenio, 
Phys. Rev. Lett. 107, 160601 (2011), arXiv:1102.1483. 

[48] M.P. Woods and M.B. Plenio, arXiv:1508.07354 (2015), 
1508.07354. 

[49] A. Go and A. J. Millis, Phys. Rev. Lett. 114, 016402 
(2015), arXiv:1311.6819. 

[50] A. Liebsch and H. Ishida, J. Phys.: Gondens. Matter 24, 
053201 (2012). 

[51] F. A. Wolf, A. Go, 1. P. McCulloch, A. J. Millis, and 
U. Schollwock, (2015), arXiv: 1507.08650. 

[52] F. A. Wolf, 1. P. McCulloch, and U. Schollwock, Phys. 
Rev. B 90, 235131 (2014), arXiv: 1410.3342. 

[53] J. A. Justiniano, “Computational aspects of three term 
recursions”. Bachelor’s Thesis, LMU Munich (2015). 

[54] M. Abramowitz and 1. A. Stegun, eds., Handbook of 
mathematical functions (New York, Dover, 1965). 


